The Free Energy Surface of Supercooled Water 
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We present a detailed analysis of the free energy surface of a well characterized rigid model for 
water in supercooled states. We propose a functional form for the liquid free energy, supported by 
recent theoretical predictions [Y. Rosenfeld and P. Tarazona, Mol. Phys. 95, 141 (1998)], and use 
it to locate the position of a liquid-liquid critical point at Tq> = 130 ± 5 K, Pqi — 290 ± 30MPa, 
and pc> = 1-10 ± 0.03 g/cm 3 . The observation of the critical point strengthens the possibility that 
SPC/E water may undergo a liquid-liquid phase transition. Finally, we discuss the possibility that 
the approach to the liquid-liquid critical point could be pre-empted by the glass transition. 
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I. INTRODUCTION 

The thermodynamic description of supercooled water 
has been a major research topic in recent years. Strik- 
ing anomalies — such as the existence of a minimum in 
the isothermal compressibility Kt along isobars, the in- 
crease of the isobaric specific heat Cp on cooling, and the 
temperature of maximum density Tm d along isobars — 
characterize the behavior of liquid water Jl]-[|. In par- 
ticular, the study of supercooled states of water sheds 
light on the understanding of the anomalous behavior of 
liquid water. The increase of Kt and C p on supercooling 
reinforces the possibility that the thermodynamic prop- 
erties of supercooled water could be different from those 
of simple liquids. Speedy and Angell proposed a scenario 
in which the increase of Kt and C p is related to a re- 
entrant spinodal line in the phase diagram of water by 
postulating the existence of an ultimate limit of stability 
for the liquid phase on cooling [Q. 

More recently, increased computing power has made 
possible the numerical study of the thermodynamic prop- 
erties of models for water. In particular, supercooled 
states, where relaxation times increase by several orders 
of magnitude over typical liquid values, have become 
computationally accessible. It has been shown that ex- 
plicit atom models (such as ST2 |L TIP4P/TIP5P [S, 
and SPC/E 0), as well as lattice ||] and continuum ||] 
models, are able to reproduce the anomalous thermody- 
namic properties of water. In all the atomistic models 
that have been studied, it has been found that the spin- 
odal line is not re-entrant. Additionally, for the ST2 
model, the existence of a novel liquid-liquid critical point 
has been directly observed in molecular dynamics simu- 
lations PQ| . Hence, it has been proposed that the anoma- 
lous thermodynamic properties of liquid water could be 
related to a liquid-liquid phase transition. According to 



this hypothesis, two distinct forms of liquid water, sepa- 
rated by a first-order transition, may exist below a crit- 
ical temperature T c ; such a critical point would account 
for the unusual increases in the thermodynamic response 
functions on cooling. Unfortunately, in water, the es- 
timated T c is below the homogeneous nucleation tem- 
perature, i.e., inside the so-called "no-man's land". This 
notwithstanding, recent experiments have probed the 
possible thermodynamic scenarios which characterize liq- 
uid water . 

From a simulation point of view, the ST2 model is the 
only one that allows a direct study of the liquid-liquid 
critical point; an increase of many orders of magnitude 
in computing power is needed for a direct detection of a 
critical point in other point charge models. Also, in su- 
percooled states at the same T and P, ST2 molecules are 
more mobile compared to real water. This feature has 
been exploited for equilibrating configurations at rela- 
tively low T jl0],[n| . The ST2 potential is over-structured 
compared to water, and the equation of state is shifted 
to higher values of pressure P, and temperature T pp| . 

Among the molecular potentials which have been stud- 
ied in detail, a significant role has been played by the 
extended simple point charge (SPC/E) model, both be- 
cause of its simplicity and its success in capturing the 
properties of water in the bulk state JlJ-jlTj], as well as 
in biological systems ]18|] . The SPC/E model has three 
point charges, located at the atomic centers of the water 
molecule. SPC/E is under structured, with its equation of 
states shifted to lower values of P and T compared to wa- 
ter [Q. Also, in the supercooled regime, at the same T 
and P, SPC/E molecules are less mobile than real water 
molecules [ 15]- 17 . Since it has been shown that the ST2 
and SPC/E models bracket the thermodynamic behavior 
of water in the T — P plane |l4|] , it would be encouraging 
to clearly detect the presence of a liquid-liquid critical 
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point also in the SPC/E potential. Unfortunately, the 
reduced diffusivity of SPC/E compared to ST2 makes it 
impossible to study directly the low T and high P region, 
where the SPC/E second critical point should be located. 

Here we propose a functional form for the free energy 
surface of the SPC/E model in the low temperature re- 
gion. Our work is supported by recent theoretical pre- 
dictions for the T dependence of the potential energy 
in supercooled states [fl9|, which have been tested for 
several model liquids |20| , |21| , p3| , p^| . The calculated func- 
tional form provides a good description of the thermo- 
dynamic quantities in the region where simulations are 
feasible, and predicts the presence of a liquid-liquid crit- 
ical point C at T c < = 130 ± 5 K, P c > = 290 ± 30MPa, 
and pc" = 1.10 ± 0.03 g/cm 3 , in reasonable agreement 
with prior estimates based on the characteristic shift 
in thermodynamics properties between the SPC/E and 
the ST2 model. 



II. THE SPC/E HELMHOLTZ FREE ENERGY 

The numerical data used to calculate the Helmholtz 
free energy F = E — TS are obtained from the long 
molecular dynamics simulations of Ref. [Gil for 42 dif- 
ferent state-points, comprising 7 different densities and 
6 different temperatures. The simulation results for the 
total energy E and pressure P are used here to recon- 
struct F in the region T > 210K, as we describe below. 
As noted in Ref. jwj, the energy £ as a function of V 
develops an increasingly pronounced convexity on lower- 
ing T. This signals the possibility of a phase transition, 
as F = E — TS will be then also convex at low T. 

Simulations of the SPC/E model below T « 200K are 
not feasible at the present time, as the time needed to ob- 
serve equilibrium metastable properties exceeds currently 
available resources. Here, the simulation data for SPC/E 
water are limited to the region T > 200K. To investigate 
the phase behavior at lower T, we exploit the recently- 
proposed relationship for the low-T dependence of the po- 
tential energy U along isochoric paths j^J. Specifically, 
the low-T behavior of the potential energy is predicted 
to follow the functional form 11911 



Uq(V), the clear negative concavity of Uq(V) at large vol- 
umes indicates that if the T 3 / 5 law would hold down to 
T = OK, then the extrapolated liquid free energy would 
imply a two-phase coexistence at zero temperature. As 
will be discussed in more detail later, at T = OK the two 
phases are separated by a first order transition around 
P = 380 MPa. 

Since the V dependence of Uo and a is smooth, we 
derive a functional form Uat(V,T), by fitting the val- 
ues of Uq(V) and a(V) with sixth degree polynomials 

Uo(V) = TLoKV" and a(V) = En=o a »^" We thus 
obtain 



Ust(T, V) = Y1 a ^ yn + T3/5 E Kyn 



(2) 



n=0 



n=0 



The a n and b n values are reported in Table BJ. We 
find almost identical values of Uq and a if we trun- 
cate Eq. (0) at order V 5 . The resulting expression 
E(V, T) = 3k B T+U m (V, T) for the total energy describes 
very well the simulation results, as shown in Fig. 

We obtain the entropy S using the thermodynamic re- 
lation 



S(T,V) = S(T ,Vo) + - 



T,V 



(dE + PdV) (3) 



T ,Vo 



where the state point (To, Vq) is some reference state 
point. We calculate the temperature-dependence of S 
along isochores from Eq. (|^), by performing thermody- 
namic integration along constant V paths. S(T, V) is 
given by 



T dT ( dE" 

K 

T\ 3 



To T \dTj v 



S{T, V) = S(T ,V) + 

= S(T Q , V) + 3k B In - ~a{V) (t"I - ) 

The unknown S(Tq,V) function can be evaluated, at 
any chosen To, from the knowledge of the ^-dependence 
of P using Eq. ||. To calculate S(T , V) we fit the simu- 
lation data for P(T — 300K, V) again with a sixth-order 
polynomial 



(4) 



U tit {T,V) = U Q (V)+a(V)T 3 / 5 . 



(1) 



Here Uo represents the T = K value of Ufa , which for a 
classical system may also be identified with F(0, V). The 
functional form of Eq. (Q) has been shown to describe 
the temperature dependence of the potential energy in 
several different models, ranging from Lennard-Jones to 
Yukawa potentials |2(i|-|24|l . Although no specific predic- 
tion has been presented until now for molecular systems, 
we find that in the temperature range between T = 200 
and T = 300 K the SPC/E potential energy is very well 
described by the T 3 / 5 law, as shown in Fig. [v|. The 
volume dependence of Uq(V) and a(V) are reported in 
Fig. and in Table | Since F(T = 0, V) coincides with 



P fit (T = 300K, V) 



ra=0 



l v n . 



(5) 



The values of the resulting c„ coefficients are reported in 
Table [n]. From elementary calculus, 



S(T ,V) = S(T ,V ) 



E(T ,V ) - E(T ,V) 
To 



(6) 



JL r T/" +1 — V n+1 



n + 1 



To 



The only unknown quantity left is F(Tq,Vq), which 
can be calculated, if needed, starting from a known ref- 
erence point (as for example an ideal gas state, as done 
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in ref. pH ) and performing thermodynamic integration 
up to Vq,T . The resulting expression for F(T,V) = 
E(T, V) - TS{T 7 V) can then be used to calculate ther- 
modynamic properties of SPC/E water. 



III. RESULTS 

First, we compare in Figs. |y| and ^ the values of 
E &t = (1 - Td/dT)F m and P fit = -dF m /dV with the 
simulation results for T < 300K. Note that the deriva- 
tives eliminate the unknown constant F(Vq, Tq). We also 
calculate the line of density maxima, Tmd, defined as 
the locus where (dV/dT)p = 0. The predicted line is 
compared with the results of the simulations in Fig. [v|. 

We next use the expression for F to calculate the ther- 
modynamic properties for T < 200K where simulations 
are not feasible. The free energy expression proposed 
depends primarily on the assumption of the T 3 / 5 depen- 
dence of the potential energy in supercooled states. The 
theoretical prediction and the quality of the T 3 / 5 descrip- 
tion reported in Fig. ^ suggests that we may meaning- 
fully extrapolate the calculation to a temperature lower 
than the one for which equilibration is feasible at the 
present time, and search for the possibility of a liquid- 
liquid critical point. 

We calculate Pfit and find that, at temperatures lower 
than 130 ± 5K, a van der Waals loop (Fig. [v|) devel- 
ops, signaling a first-order transition between two liquid 
phases. The common tangent construction |2j| for the 
Helmholtz free energy Ffi t allows us to calculate the co- 
existence line; further, we calculate the spinodal lines 
(dPftt/ 9V)t — 0. The coexistence line meets the spin- 
odal at a critical point C, which we find at Tc> = 130 ± 
5K, P c > = 290 ± 30 MPa, and p c < = 1.10 ± 0.03 g/cm 3 . 

The resulting SPC/E phase diagram is shown in Fig. ^ 
in both the (P — T) and (p — T) planes. Fig. [v|also shows 
the recently-calculated Kauzmann temperature Tk lo- 
cus defined as the temperature at which the con- 
figurational entropy vanishes |2(| . The evaluation of the 
Kauzmann locus is also based on the assumption that 
the system potential energy has a T 3 / 5 -dependence, and 
hence is fully consistent with the present free energy cal- 
culations. We note that the predicted critical tempera- 
ture is 10K below the Kauzmann temperature where 
SPC/E water is predicted to have a vanishing diffusiv- 
ity 

As a final consideration, we discuss the interplay be- 
tween the location of the critical point and the Kauz- 
mann line. Since at the Kauzmann line the configura- 
tional entropy vanishes, all equilibrium thermodynamic 
calculations lose meaning below this line. In this sense, 
the critical point in the SPC/E phase diagram should 
not be considered. In the potential energy landscape 
paradigm p6|-p9]| , the system would be trapped in a sin- 
gle basin reached at Tk- None-the-less, the free energy 
below Tk can still be calculated by considering its sep- 



arate parts. The contribution to the free energy due to 
the multiplicity of basins sampled would be fixed at its 
value at Tk, i.e. zero. Thus, below Tk the intra-basin 
free energy coincides with F. At low T, frequently, a 
model based on a harmonic solid is appropriate for such 
a calculation 20-2^]. The 'free energy calculated will still 
display a critical point (but slightly shifted compared to 
the present estimate) since the basic mechanism which 
gives rise to the low-T instability is the shape of E(V, T), 
which is already convex well above Tk- 



IV. CONCLUSIONS 

In this article, we have presented a technique of eval- 
uating thermodynamic quantities in the supercooled re- 
gion, in a T-range where equilibrium simulations are not 
feasible due to extremely long equilibration times. The 
relevant result of this analysis, applied to the SPC/E po- 
tential, is a clear indication that the free energy surface 
develops a region of negative curvature on cooling. A 
liquid-liquid critical point develops, in analogy with the 
behavior of the ST2 model, for which the location of the 
critical point is within the region where equilibrated con- 
figurations can be calculated. 

The predictions reported in this manuscript are based 
on a functional form for the liquid free energy, supported 
by recent theoretical predictions Jl^]. Of course, changes 
in the temperature dependence related to novel phenom- 
ena which may take place outside the range where data 
are available may break the validity of the extrapolation. 
In the case of real water, for example, it has been argued 
that a change in the T-dependence of the thermodynamic 
properties takes place in the "no-man's land" |50|. In the 
case of SPC/E water, if such change takes place, it must 
occur at T < 200i^, i.e. in the region where simulations 
are not feasible. This would effect our estimate of the 
location of the critical point. However, the existence of 
a region of negative curvature already in the T-region 
where simulations are feasible supports the possibility 
that the liquid-liquid critical transition would take place 
at lower temperatures, independently from the assumed 
T 3 / 5 law. 

Our results have a particular relevance, since, as pre- 
viously noticed, ST2 and SPC/E typically bracket the 
thermodynamic behavior of the real liquid. The evi- 
dence presented here that the SPC/E potential displays 
a critical point at low T and high P strengthens the 
possibility that, below the homogeneous nucleation tem- 
perature, water may undergo a liquid-liquid (or glass- 
glass) phase transition; the two distinct liquid phases 
that would appear below C could correspond to the two 
observed amorphous forms of solid water, low density 
amorphous and high density amorphous ice. Indeed, such 
a transition could be observed in the glassy state even if 
T c > <T K ,as we find for the SPC/E model. 

The thermodynamic analysis presented here also al- 
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lows us to grasp the origin of the presence of the critical 
point. Indeed, the presence of the critical point arises 
from the negative concavity of E(T, V), which for T > Tq 
is compensated by the —TS(T, V) contribution. Note 
that, as previously observed J10[ , the negative concavity 
of E(T, V) already appears in the T-region where equi- 
librium simulations are feasible, suggesting an inevitable 
phase transition as the product TS becomes progressively 
smaller with decreasing T. Such negative concavity of 
E(V,T) is also found in supercooled water pl|. 
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FIG. 1. Fit of the potential energy along isochores with 
the functional form Uq + aT 3//5 . Symbols denote different 
molar volumes. For sake of clarity, the different isochores 
are shifted by -1 kj/mol each in order to avoid overlaps. 



FIG. 3. Comparison between the energy E calculated from 
simulations |l4| and from the free energy approach described 
here. 
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FIG. 2. The volume dependence of Uo(V) from eq. (g). 
Note that this coincides with F(T = K, V). The negative 
curvature implies the presence of an unstable region in the 
phase diagram at low T. 
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FIG. 4. Comparison between pressures as calculated from 
Fig. § 
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FIG. 5. Pressure at T = 100K as calculated. The equilib- 
rium pressure is obtained by the Maxwell construction. 
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FIG. 6. P-T (upper panel) and p-T (lower panel) phase 
diagrams. The coexistence line, the spinodals and the Tmd 
line from our free energy. Squares are Tmd points obtained 
from simulations; the triangles are the Kauzmann bound- 
ary [^6| for SPC/E water below which diffusivity is pre- 
dicted to vanish p7| . 



TABLE I. Fitting parameters for the potential energy to 
Eq. (El). The fits are made for 210 K < T < 300 K. 



V (cm 3 /mol) 


Uo (kJ/mol) 


a (kJ/(mol-K 3/5 )) 


18.96421 


-83.41894 


1.1970260 


18.01600 


-80.86653 


1.1000960 


17.15810 


-78.47431 


1.0130790 


16.37818 


-76.65199 


0.9468765 


15.01333 


-74.77946 


0.8767148 


13.85846 


-74.50920 


0.8653371 


12.86857 


-74.91184 


0.8835562 



TABLE II. Polynomial fitting coefficients for Uo(V), a(V) (see Eq. 
dimensions of the coefficients depend on the term n of the expansion. 



@) and for P{T = 300K, V) (Eq. |). Note that the 



n 


a n (kJ-mor-Vcm 3 ") 


b n (kJ-mol n -7(K 3/5 -cm 3n ) 


c„ (MPa-mor/cm 3 ™) 





76.617 


-1.8261 


6.8671 x 10 4 


1 


-30.435 


0.61927 


2.4466 x 10 3 


2 


1.279.8 


-2.9301 x 10 -2 


2.8096 x 10 3 


3 


5.0719 x 10~ 2 


-9.6397 x 10" 4 


3.0755 x 10 2 


4 


4.4964 x 10" 5 


-2.63884 x 10 -5 


1.2970 x 10 1 


5 


-3.7530 x 10" 4 


1.07393 x 10" 5 


0.1871 


6 


1.1301 x 10 -5 


-3.1252 x 10 -7 


3.4974 x 10 -4 
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